# libraries
library(dplyr)
library(texreg)
library(stats)
library(MASS)
library(ROCR)
library(survival)
library(gdata)

# remove everything 
remove(list = ls())

# Data prep ####

# import the data
d <- read.csv("data.csv")

# Data: Select Variables ####
d <- d %>% dplyr::select(SideA, SideB, DyadId, GWNoA, Year, 
                         Negotiation, demref, polyarchy, duration_dy_ln, imp.gdppcln, imp.popln,
                         Intensity,
                         Region, RebStrBin.imp, Ethnic.incomp, Incompatibility, 
                         medprev2y,
                         ExChange, oda.pc.ln, UNPKO.total.ln, BdBest_ln, tsln, tsln2, tsln3)

## factor variables ####

# Incomp
d$Incompatibility <- as.factor(d$Incompatibility)
d$Intensity <- as.factor(d$Intensity)
d$Region <- as.factor(d$Region)


# Preliminaries ####
## Colors ####
percent <- 60
fark <- 60
col.neg <- rgb(red = col2rgb("gray40")[1], 
               green = col2rgb("gray40")[1],
               blue = col2rgb("gray40")[1], max = 255,
               alpha = (100-percent)*255/100)
col.both <- rgb(red = col2rgb("gray60")[1] - fark*2, 
                green = col2rgb("gray60")[1] - fark,
                blue = col2rgb("gray60")[1], max = 255,
                alpha = (100-percent)*255/100)

col.ref <- rgb(red = col2rgb("gray80")[1] - fark, 
               green = col2rgb("gray80")[1] - fark/2,
               blue = col2rgb("gray80")[1], max = 255,
               alpha = (100-percent)*255/100)

col.white <- rgb(red = col2rgb("white")[1], 
                 green = col2rgb("white")[1],
                 blue = col2rgb("white")[1], max = 255,
                 alpha = (100-0)*255/100)
col2rgb("indianred")

ct <- function(x){x*255/100}
mycol <- "lightsalmon1"
col2rgb(mycol)
mycolline <- "lightcoral"

### grey colors ####
mygrey  <- rgb(t(col2rgb("grey80")), max = 255, alpha = 125)
mygrey1 <- rgb(t(col2rgb("grey60")), max = 255, alpha = 125)
mygrey2 <- rgb(t(col2rgb("grey70")), max = 255, alpha = 125)
mygrey3 <- rgb(t(col2rgb("grey80")), max = 255, alpha = 125)
mygrey4 <- rgb(t(col2rgb("grey90")), max = 255, alpha = 125)

## mode finder function ####
mode.proper <- function(x){
  x <- na.omit(x)
  md<- as.numeric(names(table(x))[table(x)==max(table(x))])
  return(md)}

# Descriptive statistics ####

# Descriptive statistics reported in main text and Online Appendix A Table III
table(d$Negotiation, d$demref) # Appendix A Table III
apply(table(d$Negotiation, d$demref), 2, sum) # Appendix A Table III
# corresponding percentages
round((table(d$Negotiation, d$demref)[,1] / sum(table(d$Negotiation, d$demref)[,1])) * 100,2)
round((table(d$Negotiation, d$demref)[,2] / sum(table(d$Negotiation, d$demref)[,2])) * 100,2)

# Main Models Dyad-Level ####

## Table 1: Models 1-4 ####
# model 1 
fm1 <- as.formula (Negotiation ~ polyarchy + demref + duration_dy_ln + imp.gdppcln + imp.popln + Intensity +
                     Region + tsln + tsln2 + tsln3)
#model 2
fm2 <- as.formula (Negotiation ~ polyarchy + demref + duration_dy_ln + imp.gdppcln + imp.popln + Intensity +
                     Region + RebStrBin.imp + Ethnic.incomp + ExChange + oda.pc.ln + medprev2y +
                     tsln + tsln2 + tsln3 )

#model 3
fm3 <- as.formula (Negotiation ~ polyarchy + demref + duration_dy_ln + imp.gdppcln + imp.popln + Intensity +
                     Region + RebStrBin.imp + Ethnic.incomp + ExChange + oda.pc.ln + medprev2y +
                     tsln + tsln2 + tsln3 + BdBest_ln + UNPKO.total.ln)

# Main model regressions
m1 <- glm(formula = fm1, data = d, family = binomial(link = "logit"))
m2 <- glm(formula = fm2, data = d, family = binomial(link = "logit"))
m3 <- glm(formula = fm3, data = d, family = binomial(link = "logit"))

#model 4 ## fixed effects (conditional logit)
library(survival)
m4 <- clogit(Negotiation ~ polyarchy + demref + duration_dy_ln + imp.gdppcln + imp.popln +
               Intensity + ExChange + oda.pc.ln + medprev2y +
               tsln + tsln2 + tsln3 + BdBest_ln + UNPKO.total.ln + strata(DyadId), data = d)

screenreg(list(m1, m2, m3, m4))

# Quantities of Interest ####

## Model 3 Simulations####
# Coef draws
library(lmtest)
library(sandwich)
coef.draws <- mvrnorm(n = 10000, mu = m3$coefficients, Sigma = vcovHC(m3, "HC1"))


# Figure 3: FIRST DIFFERENCES #####

## polyarchy ####
poly.low  <- quantile(d$polyarchy, probs =  0.25, na.rm = T)
poly.high <- quantile(d$polyarchy, probs =  0.75, na.rm = T)

polyH <- model.matrix(m3)
var <- "polyarchy"
polyH[, var] <- as.numeric(poly.high)
polyL <- polyH
polyL[, var] <- as.numeric(poly.low)

ystar.1   <- t(coef.draws %*% t(as.matrix(polyL)))
ystar.2   <- t(coef.draws %*% t(as.matrix(polyH)))

yhat.sim1 <- exp(ystar.1) / (exp(ystar.1) + 1)
dim(yhat.sim1)
yhat.sim2 <- exp(ystar.2) / (exp(ystar.2) + 1)

fd.poly <- yhat.sim2 - yhat.sim1
dim(fd.poly)

fd.poly.pe <- quantile(apply(fd.poly, 2, mean), probs =  0.50)
fd.poly.lb <- quantile(apply(fd.poly, 2, mean), probs =  0.025)
fd.poly.ub <- quantile(apply(fd.poly, 2, mean), probs =  0.975)

## Executive change ####
EC.1 <- model.matrix(m3)
var <- "ExChange"
EC.1[, var] <- 1
EC.0 <- EC.1
EC.0[, var] <- 0

ystar.1   <- t(coef.draws %*% t(as.matrix(EC.0)))
ystar.2   <- t(coef.draws %*% t(as.matrix(EC.1)))

yhat.sim1 <- exp(ystar.1) / (exp(ystar.1) + 1)
dim(yhat.sim1)
yhat.sim2 <- exp(ystar.2) / (exp(ystar.2) + 1)

fd.EC <- yhat.sim2 - yhat.sim1
dim(fd.poly)

fd.EC.pe <- quantile(apply(fd.EC, 2, mean), probs =  0.50)
fd.EC.lb <- quantile(apply(fd.EC, 2, mean), probs =  0.025)
fd.EC.ub <- quantile(apply(fd.EC, 2, mean), probs =  0.975)


## Dem Reform #####
D.ref <- model.matrix(m3)
var <- "demref"
D.ref[, var] <- 1
D.noref <- D.ref
D.noref[, var] <- 0

ystar.1   <- t(coef.draws %*% t(as.matrix(D.noref)))
ystar.2   <- t(coef.draws %*% t(as.matrix(D.ref)))

yhat.sim1 <- exp(ystar.1) / (exp(ystar.1) + 1)
dim(yhat.sim1)
yhat.sim2 <- exp(ystar.2) / (exp(ystar.2) + 1)

fd.Dref <- yhat.sim2 - yhat.sim1
dim(fd.Dref)

fd.Dref.pe <- quantile(apply(fd.Dref, 2, mean), probs =  0.50)
fd.Dref.lb <- quantile(apply(fd.Dref, 2, mean), probs =  0.025)
fd.Dref.ub <- quantile(apply(fd.Dref, 2, mean), probs =  0.975)

## Rebel Strength ####
RebStr0 <- model.matrix(m3)

var <- "RebStrBin.imp"
RebStr0[, var] <- 0
RebStr1 <- RebStr0
RebStr1[, var] <- 1

ystar.1   <- t(coef.draws %*% t(as.matrix(RebStr0)))
ystar.2   <- t(coef.draws %*% t(as.matrix(RebStr1)))

yhat.sim1 <- exp(ystar.1) / (exp(ystar.1) + 1)
yhat.sim2 <- exp(ystar.2) / (exp(ystar.2) + 1)

fd.RebStr <- yhat.sim2 - yhat.sim1
dim(fd.RebStr)

fd.RebStr.pe <- quantile(apply(fd.RebStr, 2, mean), probs =  0.50)
fd.RebStr.lb <- quantile(apply(fd.RebStr, 2, mean), probs =  0.025)
fd.RebStr.ub <- quantile(apply(fd.RebStr, 2, mean), probs =  0.975)

## Mediation ####
Med0 <- model.matrix(m3)

var <- "medprev2y"
Med0[, var] <- 0
Med1 <- Med0
Med1[, var] <- 1

ystar.1   <- t(coef.draws %*% t(as.matrix(Med0)))
ystar.2   <- t(coef.draws %*% t(as.matrix(Med1)))

yhat.sim1 <- exp(ystar.1) / (exp(ystar.1) + 1)
yhat.sim2 <- exp(ystar.2) / (exp(ystar.2) + 1)

fd.Med <- yhat.sim2 - yhat.sim1
dim(fd.Med)

fd.Med.pe <- quantile(apply(fd.Med, 2, mean), probs =  0.50)
fd.Med.lb <- quantile(apply(fd.Med, 2, mean), probs =  0.025)
fd.Med.ub <- quantile(apply(fd.Med, 2, mean), probs =  0.975)

## Battle Deaths #####
v.BRDlow  <- quantile(d$BdBest_ln, probs =  0.25, na.rm = T)
exp(v.BRDlow)
v.BRDhigh <- quantile(d$BdBest_ln, probs =  0.75, na.rm = T)
exp(v.BRDhigh)


BRDlow <- model.matrix(m3)
BRDlow[, "BdBest_ln"] <- v.BRDlow
BRDhigh <- model.matrix(m3)
BRDhigh[, "BdBest_ln"] <- v.BRDhigh

ystar.1   <- t(coef.draws %*% t(as.matrix(BRDlow)))
ystar.2   <- t(coef.draws %*% t(as.matrix(BRDhigh)))

yhat.sim1 <- exp(ystar.1) / (exp(ystar.1) + 1)
yhat.sim2 <- exp(ystar.2) / (exp(ystar.2) + 1)

fd.BRD <- yhat.sim2 - yhat.sim1
dim(fd.BRD)

fd.BRD.pe <- quantile(apply(fd.BRD, 2, mean), probs =  0.50)
fd.BRD.lb <- quantile(apply(fd.BRD, 2, mean), probs =  0.025)
fd.BRD.ub <- quantile(apply(fd.BRD, 2, mean), probs =  0.975)

## Duration ####
v.DURlow  <- quantile(d$duration_dy_ln, probs =  0.25)
exp(v.DURlow)
v.DURhigh <- quantile(d$duration_dy_ln, probs =  0.75)
exp(v.DURhigh)

DURlow <- model.matrix(m3)
DURlow[, "duration_dy_ln"] <- v.DURlow
DURhigh <- model.matrix(m3)
DURhigh[, "duration_dy_ln"] <- v.DURhigh

ystar.1   <- t(coef.draws %*% t(as.matrix(DURlow)))
ystar.2   <- t(coef.draws %*% t(as.matrix(DURhigh)))

yhat.sim1 <- exp(ystar.1) / (exp(ystar.1) + 1)
yhat.sim2 <- exp(ystar.2) / (exp(ystar.2) + 1)

fd.DUR <- yhat.sim2 - yhat.sim1
dim(fd.DUR)

fd.DUR.pe <- quantile(apply(fd.DUR, 2, mean), probs =  0.50)
fd.DUR.lb <- quantile(apply(fd.DUR, 2, mean), probs =  0.025)
fd.DUR.ub <- quantile(apply(fd.DUR, 2, mean), probs =  0.975)

qi <- data.frame(Dref = c(fd.Dref.pe, fd.Dref.lb, fd.Dref.ub),
                 RebStr = c(fd.RebStr.pe, fd.RebStr.lb, fd.RebStr.ub),
                 Med = c(fd.Med.pe, fd.Med.lb, fd.Med.ub),
                 Dur = c(fd.DUR.pe, fd.DUR.lb, fd.DUR.ub),
                 Brd = c(fd.BRD.pe, fd.BRD.lb, fd.BRD.ub))



dev.off()
png("First_Dif.png", units="cm", 
     width = 16, height = 10, res = 300)
par(mar = c(1,5,2,1), mfrow = c(1,1), mgp = c(3,0.8,0))
plot(seq(0.5,5, length.out = 100), seq(-0.05,0.40,length.out = 100), t = "n",
     axes = F, xlab = "", ylab = "Change in Pr(Negotiation)", main = "Mean First Differences (Model 3)")
abline(h = seq(0, 0.40, 0.05), col = mygrey3)
# Democratic reform
segments(y0= qi$Dref[2], y1 = qi$Dref[3], x0 = 1, x1= 1, lwd = 2)
points(1, qi$Dref[1], pch = 16)
# Rebel Strength
segments(y0= qi$RebStr[2], y1 = qi$RebStr[3], x0 = 2, x1= 2, lwd = 2)
points(2, qi$RebStr[1], pch = 16)
# Mediation
segments(y0= qi$Med[2], y1 = qi$Med[3], x0 = 3, x1= 3, lwd = 2)
points(3, qi$Med[1], pch = 16)
# Duration
segments(y0= qi$Dur[2], y1 = qi$Dur[3], x0 = 4, x1= 4, lwd = 2)
points(4, qi$Dur[1], pch = 16)
# BRD
segments(y0= qi$Brd[2], y1 = qi$Brd[3], x0 = 5, x1= 5, lwd = 2)
points(5, qi$Brd[1], pch = 16)
axis(2, las = 1, at = c(0,0.10,0.20,0.30,0.40),
     labels = c("0.00", "0.10", "0.20", "0.30", "0.40"), 
     tck = -0.02)
axis(2, las = 1, at = c(0.05,0.15,0.25,0.35),
     labels = NA, 
     tck = -0.01, lwd = 0.5)
axis(2, las = 1, at = seq(0,0.40, 0.05), labels = NA, tck = 1, lwd = 0.1, col = mygrey)
cxt <- 0.8
ln1 <- -1.5
mtext("Democratic", side = 1, line = ln1 , at = 1, cex = cxt)
mtext("Reform (t-1)", side = 1, line = ln1 + 1, at = 1, cex = cxt)
mtext("Rebel", side = 1, line = ln1 , at = 2, cex = cxt)
mtext("Strength", side = 1, line = ln1 + 1, at = 2, cex = cxt)
mtext("Prev.", side = 1, line = ln1 , at = 3, cex = cxt)
mtext("Mediation", side = 1, line = ln1 + 1, at = 3, cex = cxt)
mtext("Duration", side = 1, line = ln1, at = 4, cex = cxt)
mtext("Battle", side = 1, line = ln1, at = 5, cex = cxt)
mtext("Deaths", side = 1, line = ln1 + 1, at = 5, cex = cxt)
dev.off()